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ABSTRACT 


The appleeatton of analysis lattice filters to the 


problem of determining the input to a—esystem from 
observations of the system’s output (i.e., deconvolution) is 
discussed. Both linear and nonlinear systems are 
considered. Lattice filter modeling algorithms (Levinson 


and Schur) are presented. 

The theory of least-squares inverse filters is reviewed. 
This leads to a discussion of the lattice filter, which in 
ern leads to the Generalized - Lattice Theory. The 
Generalized Lattice Theory is then used to develop a 
nonlinear lattice structure. Simulations show that the 
nonlinear lattice is an effective inverse filter for both 


linear and nonlinear systems. 
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I. INTRODUCTION 


The problem of estimating a signal based upon 
observations of a related signal is one of the most 
important operations in signal processing. The input and 
output Signals of a linear system are related by the 
convolution operation 

t 
y(t) = h(t) * x(t) | h(t-?T)x(T)d¥ eee 
— 0 


where h(t) is a causal, linear time-invariant (LTI), system 


impulse response. Since this thesis deals primarily with 
discrete digital signals, continuous time signals are 
sampled at uniform intervals, T, and are represented by 
discrete sequences x(nT) = x(n) for n = O,1,2,...,N. 


Discrete convolution for a LTI system is defined by 
n 
y(n) = h(n) * x(n) = )) h(n-m)x(m) (1.2) 
m=-@ 

and is shown in Figure 1.1. The basic deconvolution 
problem is to estimate the signal x(n), assuming that both 
y(n) and h(n) are known. Figure 1.2 depicts the inverse 
filtering process of recovering the input signal from the 
output signal by removing the system’s impulse response. 

Deconvolution has important applications ina variety 


of fields: Radar, communications, image processing, speech 


Original Observed 
Signal Signal 
SYSTEM 


x(n) y(n) 


h(n) 


y(n) = x(n) * h(n) 


Figure 1.1: Convolution 


Observed Original 
Signal Signal 






INVERSE FILTER 


y(n) x(n) 


g(n) 


x(n) = y(n) * g(n) 


Figure 1.2 Deconvolution 


synthesis, seismology. For example, in image processing, 
deconvolution is used to recover the representation of the 
original object, x(m,n), from the representation of its 
image, y(m,n), by removing the blurring caused by the opti- 
cal system’s point-spread function, h(m,n). A common 
problem which arises in geophysics involves the deconvolu- 


tion of a seismic trace, y(n), into the approximately known 


impulsive waveform, het) and the desired reflection 
response, x(n), which reveals the structure of the layered 
Earth. In other applications, h(n) may represent’ the 
impulse response of a transmission channel, magnetic 


recording medium, or measurement device which broadens’ and 
smears (intersymbol interference) the desired message x(n). 
There have been numerous approaches to the linear 
deconvolution problem, including least-squares filtering, 
linear inverse theory, linear programming, and homomorphic 
Signal processing. The first part of this thesis deals with 
least-squares filtering techniques; the application of 
Kalman, waveshaping, and lattice filters to deconvolution is 
reviewed. Next, the lattice filter discussion is extended 
to the theory of the generalized lattice filter. Finally, 
nonlinear system theory is’ briefly reviewed, and the 
nonlinear lattice filter is developed and applied to the 
inverse filtering problem. Computer simulation results for 


both the linear and nonlinear lattice filters are presented. 


II. LINEAR LEAST-SQUARES DECONVOLUTION TECHNIQUES 


A. INTRODUCTION 
This chapter begins with a discussion of the principles 


of least-squares filtering theory since deconvolution is 


essentially an inverse filtering process. Three specific 
types of least-squares filters are then introduced: Kalman, 
waveshaping, and lattice filters. The application of each 


of these three filter types to the problem of linear decon- 
volution is studied. Finally, the chapter concludes’ with 
the presentation of computer simulation results obtained for 
deconvolution experiments using the lattice filter. 

The goal of deconvolution is to recover the input 
Signal, x(n), to a system based upon observed values of the 
system’s output signal y(n). An optimal processor must be 
determined to - produce the best possible estimate of x(n) 
based upon present and past values of the output y(n). A 


traditional measure for defining the "best" signal processor 


is the minimum mean-squared error criterion. Using this 
criterion, the estimate x(n) is defined by a linear 
combination of the observed values yaén ) . Assuming 
causality, and that the observed signal is windowed to 


include only the M past samples, X(n) is written as 
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A . ° 
x(n} = ys hing.) yt 1) eee 2. |) 
1=n-M 


where h(n,i) = O for n < 1. The estimation error is given by 
ein) = x(n) = Z(n). To find the optimum signal processor, 
the h(n,i) coefficients which minimize the mean-square 


estimation error J, where 
2 7 2 

J = Efe (n)] = E{ (x(n) - x(n)) J; 2-2) 
must be determined. This is known as the Wiener filter 
formulation of the problem. [Ref. l:pp. 113,116] 

The least-squares theory of filtering began in the 
1940’s with the work of Norbert WIENER [Ref. 2:pp. 147-148]. 
WIENER developed a frequency domain procedure to design 
optimum filters, where optimality was defined by minimiZing 
a mean-square error performance criterion. The Wiener 
filter is conventionally applied to linear time-invariant 
systems with stationary statistics when it is desired to 
separate one signal from another. In the early 1950’s, the 
Wiener filter was extended to include time varying and 
nonstationary statistics, but the calculations are 
cumbersome [Ref. 3:p. 1]. 

The mean-square estimation error J({(n) 1S minimized by 
setting its partial derivatives, with respect to each of the 


filter coefficients h(n,i), equal to zero. Thus: 
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QJ (n) 


m—_ | =" we (2. on 
Ghi(n,i) 
for 1.2 n-M, nM... on This yields a set of Mtl linear 


Simultaneous equations, called the orthogonality equations 


where 

R (n,i) = Efe(n)y(i)J = 0 (2943 

ey 
for NeMOPGVes “Seon and no 9a oc. Ron | ees the 
ey 
crosscorrelation function between the error. signal, e(n), 
and the data y(n). Substituting e(n) = x(n) - Xn) into the 
orthogonality equations gives the normal , or Wiener-Hopf 
equations: 
n 
E{x(n)y(i)] = )) h(n, k)ELy(k)y(i)] (2.5a) 
k=n-M 
or 
n 

R (nya) =). eo (ic) Gk eee (2.5b) 

Sy. k=n-M yy 
Since the autocorrelation function R of the input signal 

yy 
and the crosscorrelation function R of the desired output 
XY 

Signal with the input signal are known quantities, the M#l 


equations can be solved for the optimal filter weights 
nn, 1) 4-2 fan. ee If data vectors are defined so that 


T 
x(n) = [x(n=M), x(n-Mt1)4 22 coae (2.649 
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and 


Aly 
wie | Viti) ey (n—M+1.. , yen) | 2.6 bs) 


then the M+1l Wiener-Hopf equations can be written as 


T T 
E{x(n)y (n)] = hEly(n)y (n)] (2-7) 
where h = [{hi(n,n-M), h(n-M+l1),...,h(n,n)]. Assuming that 
the signals x(n) and y(n) are stationary, then the filter 
coefficients are time-invariant and h(n,k) = h(n-k); and 


equation (2.7) can be written in matrix form as 


Ryy(QO) Ryy(1l1) Ryy(2) ... Ryy(M) h(O){-—-  4Rxv(0) 
Ryyil) Ryy(0) Ryy(1) ... h(1) Rey41) 
Ryy(2) Ryy(1) Ryy(Q) ... nz) aro y tc) (Zee oO) 
Ryy (M) Ryy(0) h(M) Rxy(M) 


Now, the optimal coefficients can be solved by inverting 
the autocorrelation matrix, or by exploiting the matrix’s 
Toeplitz structure (all elements are the same on any given 
diagonal ) to employ the more efficient Levinson algorithm. 
(Levinson’s algorithm will be discussed in the development 
of the analysis lattice filter.) 

The minimized value of the mean-square estimation error 
can now be computed, and is found to be 

Z n 


oy (n) = E[{x (n)] - % h(n,i)Ely(i)x(n) ] (2.9) 
min iz=n-M 


13 


T li -1 
= R (0)-E(x(n)y (n)JEly(n)y (n)] Ely 
xX ; 
(Ref. 4:p. 148]. These values correspond to the diagonal 


elements of the covariance matrix of the estimation error, 


T 
R = Efe(n)e (n)] (2 Ton 
ee 


A 
where e(n) = x(n) - x(n) (Ref. 1l:p. 120]. Furthermore, the 
estimate of x(n) is given by 


R T al -1 
x(n) = E{x(n)y (n) JEL¥y(n)ys ie (2.11) 


uN 


hy(n) 


This estimate can be thought of as the projection of the 
desired signal x(n) onto the space spanned by the components 
of the observation vector y(n). The minimized estimation 
error vector 1s orthogonal, or normal, to the estimate 
x(n). (Ref. 4:p. 147] 

This completes the overview of least-squares filtering. 
The following sections present several linear deconvolution 
techniques which employ this criterion: Kalman filtering, 


spiking filters, and lattice filters. 


B. KALMAN FILTER 
In the early 1960’s, R.E. KALMAN introduced an optimal 
recursive filter based on state-space time-domain methods 


[Ref. 2:pp. 267-268]. The Kalman filter estimates the state 
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of a linear system, and is optimal in the sense that it min- 
imizes the mean-square error of the state estimate. The 
Kalman filter is useful when the system is defined by state 
Space equations: The system signals are represented by 
random processes, and the data observations are 
contaminated by noise. The Kalman filter algorithm processes 
measurement data, and requires a priori state space models 
(known or assumed) of the system and measurement dynamics. 
Also, the statistics of the system input and measurement 
noises, as well as initial condition information, are 
required to produce the state estimate. This process is 
depicted in Figure 2.1. Here, the discrete Kalman filter 
equations will be presented, and then their application to 
deconvolution will be described. 

The state space representation of the discrete system 
and measurement models (see Figure 2.2) are written as [Ref. 


2:pp. 195-200): 


x(k) = F(K-1)x(kK-1) + w(k-1) (ace) 

z(k) = H(k)x(k) + v(k) (2.13) 
where 

x(k) = (n x 1) system state vector 

F(k) = (n x n) transition matrix 

w(k) = (n x 1) system noise vector 

z(k) = (m x 1) measurement vector 
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ayeUITISY 9384S = (3)x 


UOTJBAIISqGO = (3)z 


a}yBYS wiayshS = (})x 


MALTA 
LNGINGTOSVAWN WaLSAS 
(3)% NVWIVa = | (3)2 (3)x 
MoTzyeUIIOJU]T IOILY IOIIY 
TIO W JUSJUIZINsS evap U13a}ysSAS 


“ee R OTH EN PR FP RITD PAP AE LIDIA CNEN ITM MAAN OM) IIN 


Figure 2.1 System, Measurement, and Kalman Filter 
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Figure 2.2 System and Measurement Model 
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FWCErrtelLlyUVUUC LS AAI We Vv GIN a. 7X4 bo 8 Bw ee 


H(k) = (m x n) observation matrix 


v(k) = (m x 1) measurement noise vector. 


(Note that the time index k is used here to be consistent 
with the discrete Kalman filter literature. ) 

The noise vectors w(k) and v(k) are assumed to have zero 
mean, White, Gaussian distributions with covariance matrices 
of Q(k) = siuthia eo and R(k) = Eivikis fi 
respectively. Additionally, wandv are uncorrelated so 
that res = 0 for all k and j. The state estimation 


error is defined by 
Aw 
e(k{k-1) = x(k) - x(k{k-1) (2. tay 


and the associated (n x n) error covariance matrix 1s 


7 

P(kt(k-1) = Efe(kik-l)e (k{k-1)]. (2, bow 
An updated, a posteriori estimate of x(k) is obtained from 
the measurement z(k) and the a priori. state estimate 
R(k{k-1) by 

A A NN 

x(k) = x(kKt{k-1) + K(k) (z2(k)-H(kK)x(kKI/kK-1)] (2.16) 
where K(k) is the (n x m) Kalman gain matrix. The a 


posteriori error is given by 
A hoes 
e(k) = x(k) - x(k) (2.07 


and the associated a posteriori error covariance matrix is 
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T 
Pik) — Bletkje Ch) ]. Ze S) 
The mean-square estimation error criterion is minimized 


when 
T T -1 
K(k) =P(k{tk-1)H (k)(CH(k)P(kIkK-1)H (k)+ R(K)] (Zao) 
This optimal value of the Kalman gain matrix minimizes’ the 
individual terms along the main diagonal of Ptk). 
Substituting the optimal gain matrix K(k) into the 


expression for P(k) results in 


P(k) = (I - K(k)H(k)) P(ktk-1) , (2.20) 
where JI is the identity matrix. In order to compute 
equations (2.16) and (2.19) recursively, the ae priori 


estimates &(k+11k) and P(k+1l{tk) must be determined at time 


k. The a priori estimates are given by 
K(k+14k) = F(k)X(k) (2021) 
T 
P(k+1/k) = Efe(k+tlftk)e (k+1tk)] (Zeac} 


Ib 
Pi cieth tw oer i kpe(k)twik)) | 


T 
F(k)P(K)F (k) + Q(k) 


The Kalman filter algorithm is implemented by recursive- 
ly computing equations (2.16), C24. 19), CZ cOjtmoocl )y and 


2.22). Figure 2.3 is a diagram of the Kalman filter. The 


is 


(T-H)X 


(T-4| 1)X 





Figure 2.3. Discrete Kalman Filter 
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algorithm is initiated with the initial conditions 


Z(0) = E[x(0)] (2.23a) 


and, 


TT 
E{ (x(0)-x(0))(x(0)-x(0)) J. (2.23b) 


ECO) 


In the event that a controlling input or a deterministic 
disturbance u(k) is applied to the system, the only change 
in the above algorithm is the state model and the a poster- 


iori estimate. They become [Ref. 5:p. 130] 


pom P(e tyx(k-1) + wik-1) + Glk-1)u(k-1 ) (2.24) 


x(k) = F(k-1)2(k-1)+G(k-1)u(k-1) (2.25) 


+#K(k)[z(k)-H(k) F(k-1)%(k-1) ]. 


where G(k) is a (n x q) matrix and u(k) is a (q x 1) vector 
of input signals. 


The problem of estimating a desired signal based _ on 


noisy data observations pertains to- such fields as 
communications, controls, and geophysics. However, the 
Kalman filter can also be applied to these deconvolution 
problems. AS an example, the discrete Kalman inverse filter 


will now be applied to an exploration seismology problem. 
Typically, in the search for underground oil and gas 
deposits , a vibratory signal source generates a pulse of 


energy which is’ transmitted into the earth. In modern 
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seismology, ‘the shape of this source wavelet can be 
carefully controlled; it is chosen so that it contains only 
those frequencies which are transmitted best by the earth. 
As the source wavelet propagates through the earth, it 
encounters many different layers with various acoustic 
impedances. At these layers, both partial reflection and 
refraction occur, creating numerous transmission paths. The 
received seismic signal at the surface is composed of many 
overlapping reflected wavelets. Therefore, the seismic 
trace can be represented as the convolution of the original 
source wavelet with an impulse train representing the 
various layers of the earth. Moreover, the seismic trace is 
contaminated by measurement noise and by the phenomena of 
ghost reflections and reverberations. [Ref. 6:pp. 14-15] 


The seismic trace can be described mathematically by 


Zi) = oS Cea 1 ) ero) toe Ve) [sce styeenan + v(t) (23268 
t 
where z(t) = measured seismic trace 
s(t,T) = finite duration, time varying wavelet 
r(t) = reflectivity function of earth’s structure 


measurement noise. 


v(t) 


The seismologist must extract the structure of the earth, 
r(t), by analyzing the noisy seismic data z(t). This 
process of removing the wavelet shape from the trace and 


leaving behind the impulse train representing the reflected 
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wavelet’s strength and arrival time is that Ol 


deconvolution. CRUMP [Ref. 3ipp. 6-7] shows that if the 
seismic signal is sampled at discrete, uniformly spaced 
intervals, and if s(t,T) and r(t) are assumed to be causal, 


then z(t) is represented by 


[VIe 


Z(k) = [s(k,k-itl)r(k-it1)] + vik) CZrgt) 
i=1 
where the sample number k = i. 2 ee. and 
L = length of the wavelet given in number of samples 
=e tor Ties 
= L for k>oL . 


When M traces of K samples in length are available for 
processing then z(k) becomes a (M x 1) vector where the j-th 


component is given by 


Z (k) = y tH (k))r(k-itl)J] + v (k) (Zng57) 
i=l ji J 
meee |, 2e...,M and k = 1,2,.+..,K-. This assumes that the 
reflectivity function is the same for each trace, while the 
shape of the exciting wavelet may vary from trace to trace. 
The time-varying wavelet values are contained in the (M x L) 
matrix H(k) where the j-th row contains the L samples of the 
wavelet which generates the j-th trace: 

H (k) = s (k,k-itl) . (24:29) 

Ji J 


In vector form, the equations become 
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Z(k) = H(k)x(k) + v(k) (2.30) 


where the state, measurement, and noise vectors are given by 


ic 
x(k) = TrCRy src yy. oe ee (2.31a) 
T 
Z(k) = [z (k),2Z (K)j. 32, Zhe (2.30 
1 Z M 
T 
vik) = [v (k),v-(K)5 \.3ey Vie ; (2.30en 
t Z MS 
This is the Kalman filter measurement model. Now the state 


model must be determined. 

The state model is arrived at by assuming a general 
relationship for the reflection coefficients of x(k). The 
assumed relationship is the autoregressive equation 

L 
r(k) = ) [b (k-1)r(k-i)] + w(k-1) . (2.30 
lac) a 
Comparing equation (2.32) with the state vector x(k) yields 


the state model 
x(k) = F(k,k-1)x(k-1) + w(k-1) (2.300 


where the (M x L) transition matrix and the (M x 1) system 


noise vector are given by 
b (k-1) b (k-1) 
2 


» (2.34a) 
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T 


Pic lj wi k=1) ,0,0 , <s2,0'] ; (2.34b) 
I is the identity matrix, and 0 is the null vector. 
Equations (2.30) and (2.32) provide the measurement and 


system models for implementing the deconvolution via the 
Kalman filter. CRUMP discusses methods by which to obtain 
numerical values for the reflection coefficient vector b(k) 
and the time-varying wavelet sample matrix H(k) [Ref. 3:pp 
8-11]. Once these matrices are determined, the recursive 
Kalman filter removes the effects of s(t,T) from r(t) and 
generates the state estimate x(k) which provides L samples 
of the desired reflectivity function at each time k. 

It is interesting to note that both the Kalman and 
Wiener filters are minimum mean-Square error estimators, 
both require the same a priori knowledge of the process’~ to 
be estimated, and that both yield identical estimates. 
However, the Kalman filter does have distinct advantages 
over the Wiener filter. First, due to the matrix form of 
the state Space equations, the Kalman filter has 
multichannel capability and is equivalent to a bank of 
optimal estimators. Moreover, the Kalman filter is ideally 
Suited to computer implementation due to its discrete and 
recurSive characteristics. [Ref.2:pp. 268-269] 

The discrete least-squares approach which follows is a 


viable alternative to the Kalman £yviter algorithn, 
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particularly when state space model equations are not 


available or applicable. 


C. LEAST-SQUARES INVERSE FILTER 

When a linear system, H(z), is excited by the an input 
signal x(n), the output y(n) is defined by the convolution 
relationship y(n) = x(n) * h(n). Deconvolution involves 
finding the inverse filter G(z) such that H(z)G(z) =a 
(Note that the discrete time domain is related to the fre- 
quency domain through the equation z = =. » Where T 18 
the discrete sampling interval.) This condition transforms 
to h(n) * g(n) = a(n) in the discrete time domain; h(n) and 
g{(n) are the impulse responses of the filters H(z) and G(z), 
respectively, and ad(n) is the unit impulse function. If 
this condition is met, then the original input signal x(n) 
1s recovered at the output of the inverse filter. 

The transfer function of the causal system is defined by 


the infinite series obtained by taking the one-sided Zz- 


transform of the system’s impulse response: 
oO -i 
H(z) = Y(z)/X(z) = ) h(i)z. (2.35) 
1=0 
H(z) can be determined by inserting a known sequence x(n) 
into the system, measuring the output sequence y(n), and 
then manipulating their z-transforms. The inverse filter 


G(z) is then computed by carrying out the polynomial 


division 
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zea Ont) Zt ( 2 )iZe +e | (2.36) 


Cae Z)) c= 
-1 -2 -M 
ee 09) Me el )z + g(2)z fae er Miz +o... 
and truncating to Mtl terms if necessary. If the exact 


inverse filter iS approximated by truncating G(z) to order 


M, then its impulse reponse is given by the sequence g(n), 
mao 1,2,...,M. Furthermore, -if the original system’s. 
impulse response is represented by h(n), n=0,1,2,...N, then 
M 
amie s(n) * h(n) = d, g(m)h(n-M) fo) 
m= 


for 0 < n < NtM. This approximation to the impulse function 

improves as the order M of the inverse filter is increased. 
Now the stability of the inverse filter will be 

addressed. If H(z) has all its zeros inside the unit circle 


in the complex z-plane, it is referred to as a minimum-delay 


polynomial; the corresponding sequence h(n) is called a 
minimum phase-lag sequence. This is a sufficient condition 
to guarantee that H(z) has a stable inverse, because the 
zeros of H(z) become the poles of G(z), and G(z) is a stable 
filter if all its poles lie within the unit circle. 
Maximum- and mixed-delay Signals are obtained by 
transforming the zeros of H(z) from Zo to je where the 
Superscript ’*’ represents the Syne Been Oe operation. 


A maximum-delay sequence has all of its zeros outside the 


7 Fi 


unit circle, while the mixed-delay sequence has zeros inside 


and outside the unit circle. 

If H(z) has N zeros, then transforming these zeros 
results in at most - distinct sequences. Of note, each of 
these minimum, maximum, and mixed-delay signals have the 
Same magnitude spectrum: aie = enna [Ref. lip. 98]. 
However, they do have distinct phase spectra [Ref. 4:p. 


175]. The maximum-delay polynomial can be written as 


R * x -l X -N 
H (z) = h (N) + h (N-1)z +...¢+ h (O)zZ (2.338 


R 
The so called reverse polynomial H (zZ) is a conjugated, 


reflected, and shifted version of H(z). The corresponding 
R * * x 
maximum phase-lag sequence is h = {fh (N),h (N-1),...,h (0)} 
[Ref. 6:p. 72] While the minimum-delay filter has a causal, 
stable inverse consisting only of a memory function, the 
maximum-delay filter has an inverse which consists only of a 
stable, noncausal, anticipation function. The stable 
inverse of a mixed-delay function consists of both memory 
and anticipation functions. Filters with nonvanishing 
anticipation components are noncausal; they cannot work in 
real time since the future values of the filter input are 
not available for processing. This problem can be 
circumvented if the entire signal is first recorded prior to 


analysis; then the required future input data is available. 


[Ref. 4:p. 87] 


28 


The energy distribution in minimum, mixed, and maximum 
delay signals will now be examined. Since each of these 
Signals have an identical magnitude spectrum, they also 
have the same total energy. However, although the total 
energy is the same, the rate at which the energy builds up 
differs for the various sequences. Parseval’s theorem 


states that the total energy ina signal is given by 


Tw Z M Z 
[H(w) | dw/(2n) = ) Jh(m)I. (2323) 

~1 m=Q0 
If the partial energy is defined as 

n 2 
P(n) = ) Jh(m)t , (2.40) 
m=0 

then it can be shown that the energy builds up quickest in 
the minimum-delay sequence, and that it builds up the 
slowest in the maximum-delay sequence [Ref. 6:pp. 75-76]. In 


other words, the minimum-delay signal makes itS impact as 


soon as possible since itS energy 1S concentrated at the 


front of the sequence. The maximum-delay signal makes its 
major impact at a later time since its energy is 
concentrated at the end of the sequence. The energy curves 


associated with all the possible mixed-delay signals lie 
between these two extremes. Finally, it can be shown that 
the convolution of two minimum-delay sequences with one 
another results in a minimum-delay sequence. The convolution 


Of maximum-delay Signals results in a maximum-delay 


Uae) 


hn 
uh 


sequence. Convolution involving any other combination of 
sequences results in a mixed-delay sequence. Of note, the 
resulting maximum and minimum-delay sequences are the 
reverse of each other. [Ref. 6:pp. 73-74] 

The method described above of finding an approximate, 
finite length inverse filter consisted of simply truncating 
the exact inverse filter found by polynomial division. et 
was seen that the inverse filter G(z) attempted to transform 
the impulse response of H(z) into a unit impulse located at 
the origin. This can be thought of as an attempt by G(z) to 
undo the blurring effect of H(z) (i.e., H(z) “blurs” the 
impulse d(n) into the impulse response h(n)). If the input 
to H(z) is designated x(n), and.if the output of G(z) is 


R(n), then the error of the approximated inverse filter is 
e(n) = x(n) - K(n) for 0 < n ¢ N+M. (2.41) 


The error energy is defined by 
N+M 2 
Jit » e (n) (2 ae 
n=Q0 
For the polynomial division / truncation method, J decreases 
as the order M of G(z) increases [{Ref. 6:p. 136]. Seeking to 
minimize the error energy J with respect to the inverse 


filter coefficients leads to another approach for finding an 


approximate inverse filter. 
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Filters which minimize the mean-square error J are 
called least error energy or least-squares inverses. In 
general, the desired output sequence x(n) of the inverse 
filter G(z) could be of any shape. G(z) is then called a 
waveshaping filter. When applied to deconvolution, the 
desired output of the inverse filter is a unit impulse 
d(n-i), where i = 0,1,2,...,N+M defines the lag of the digi- 
tal inverse filter. In this case, G(z) is called an i-th 


delay spiking filter since it tries to condense the system 


impulse response h(n) into a spike with i delays. Since i = 
0,1,2,...,N+M, there are N+M+l1 possible spiking filters. As 
will be discussed, there are preferred values of the lag i 


for minimizing J, depending on the phase characteristics of 
h(n). The spiking filter problem is shown in Figure 2.4. 

It is convenient to restate the deconvolution problem in 
a matrix form based on the Yule-Walker, or autocorrelation, 
method [Ref. Ilipp. 243-245]. This is accomplished by 
defining the (M+N+1l x M+l1) system output matrix /Y, the 
(M+1 x 1) impulse response vector g, and the (N+M+1 x 1) 


input vector x as follows: 
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{<----- M ZEROS ------ > 
y(0) 90 0 nara 0 
lO} 0 ae 0 
mez ey CL) yO) ous 0 
; j : ee ¢ ¢ 0 
O y(N) yAN— 1) 
0 0 y(N) 
¢ ® 0 ° 
0 0 0 y(N) 
T 
meee), e2( 1), «2. ,8)  , (2435) 
T 
mae xO). x( loo... ,x(NtM) | : (254 3¢C) 
The above notation is for the general case of the 
waveshaping filter. For the special case of the spiking 
filter, . the y(n) components of the system output matrix Y 
are replaced by the impulse response h(n). Also, the 
desired signal x reduces to an impulse d(n-i). Now 
equations (2.37), (2.41), and (2.42) can be rewritten as 
= = “yo ; (2.44a) 
aA 
e =x - xX , and (2.44b) 
T 
Jz>ee (2.44c) 
As discussed in the introduction to this chapter, the 


least-squares criterion is satisfied by minimizing J with 
respect to the inverse filter coefficients G. This results 


in the normal equations 
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VeYe = “Vox (2.45) 


T 
which can be- solved for g. By recognizing that Y Y is 


equivalent to the (M+1 x M+1) sampled autocorrelation matrix 
Tr : 
R = Efyy ] where the M+1 length data vector is y = { y(0), 
T a 
y(1),..-.,y(N),0,0,...0] , and that Y xX is avlene@ thee 


cross-correlation column vector r, it can be seen that 


T- =) 5 -1 
= (lene a) er eel (2.46) 


-1 
where R can be evaluated efficiently by Levinson’s 


algorithm. The actual filter output is then 


-1 T 
& = Yg = (YR Y )x = Px (2.47) 

T -1T 
where P = Y(Y Y) Y is a square (N+M+1) dimensioned matrix 
called the projection or performance matrix. The filter’s 


performance improves as P approaches the identity matrix. 


Now, the estimation error and cost function are written as 


(2.48) 


| 
i> 
i 
I 
IH 
I 
& 


= pos 
T a 
Jsees= x (1 - P)x (2.4 oF 
Since in the deconvolution problem x is the spike with 1 
delays, the inverse filter output x is actually the i-th 
column of the P matrix. Also, Since g = oma the 
coefficients of the i-th spiking filter are contained in the 


-1 T 
i-th column of the matrix R Y . Moreover, the energy 
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error of the i-th spiking filter reduces to J = 1 = P(i,1). 


Therefore, in order to realize the best inverse function 
(i.e., minimize J), select the delay i for which P(i,i) is 
largest. For a chosen lag i, the filter’s performance also 


improves as the order M of the inverse filter is increased. 
Now let J(i) represent the estimation error for the i-th 
spiking filter. The grand sum of squared errors is’ defined 
as V = J({(0) + J(1) +...+ J(M+#N). It can be shown that V = 
(M+N+1) - (M+1) = N, where N is the order of the system H(z) 
[Ref. 4:p. 198]. Therefore, V is independent of order M. 
For sufficiently long spiking filters, the optimal value 
of the delay i depends on the phase characteristics of the 
Signal h(n) and the choice of the lag is governed by the 
following rules. If h(n) is a minimum-delay input, the 
spiking filter should have zero delay, i = 0. This _— 
that for a signal with its energy concentrated towards’ the 
front of the sequence, it is easiest to condense it to a 
unit impulse at the origin. If the signal iS a maximum- 


delay input, the maximum-delay spike i = N+M+1 should be 


selected. If h(n) is a mixed-delay signal, then the i 
corresponding to the largest P(i,i) should be chosen. [Ref 
Ome p. 152) 

Under certain conditions, the estimation error will go 


to zero as the length of the inverse filter tends to 


infinity. First, as previously discussed, if h(n) isa 
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minimum-delay sequence then an exact inverse can be found 
through polynomial division. The zero-delay spiking filter 
approaches this exact inverse filter G(z) as the number of 
terms M+l increases. Therefore, the estimation error goes 
to zero as M goes to infinity. The second condition is, if 
h(n) is not a minimum-delay sequence, J will approach zero 
as 1/M if the lag i of the spiking filter is chosen to be 


sufficiently large. [Ref.4:pp. 200-201] 


Up until now, the effects of measurement noise and 
imperfect knowledge of the distorting function H(z) have 
been neglected. If noise is introduced, then the output 


y(n) of the system H(z) driven by input signal x(n) becomes 
yin} <= htn)*x(n) “+ y(n) (2. 510m 


where v(n) iS assumed to be zero-mean white noise with 
variance Q. If this signal is passed through the previously 


determined inverse filter G(z), the filter output becomes 


R(n) = g(n)*y(n) = g(n)*h(n)*x(n) + g(n)*v(n) (2.51) 


it 


d(n)*x(n) + u(n) = x(n) + u(n) 


where u(n) = g(n)¥v(n) is the filtered noise signal and d(n) 
is the unit impulse function. The variance of u(n) is: 
2 M 2 T 
E{u (n)}] = Q@ ) g (n) = @gg (2.52) 
n=0 


{[Ref. l:p. 248]. This variance may be larger than _ the 
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Original noise variance Q. For example, if h(n) is a_ low 
frequencey signal, then g(n) must be very spiky in order to 
compress h(n) into a high frequency content spike. As a 
result, g(n) will likely have values greater than one. 
Therefore, the variance of u(n) will be larger than Q. This 
further degrades the estimate Qin). 

To compensate for this filtered noise, the minimization 


eriterion is modified so that 


N+M 2 r 
J= P} (d(n) - g(n)th(n)) + AQg g (225s) 
n=0 
where A is a positive parameter. The first term of the cost 


function tries to produce a good inverse filter whereas’ the 
second term tries to reduce the output noise. If Ais large, 


noise reduction is emphasized at the expense of obtaining a 


good inverse function. A small A emphasizes finding a good 
inverse function g(n), and there is little output noise 
reduction. 

Using this new cost function, the resulting normal 


equations are given by 
T T 
Pray etVAQGL)e = Y x (2.54) 
Comparing this with equation (2.45) reveals that the main 
diagonal elements of the sampled autocorrelation matrix R 
have been modified by an additive term: R(0) becomes R(0O) + 


AQ. If the Backus-Gilbert "prewhitening" parameter epsilon 
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iif 


Ri 


is introduced, where 
= NYP AUOS ; (2.508 


then the diagonal elements are written as (1 +€)R(0). Even 

very small values for the BG parameter have ae beneficial 
TT 

effect of stabilizing the inverse of the matrix (Y Y + AQI). 


{[Ref. lip. 246] 


Da SLATTECE “FITTER 
i TS INncroduetlon 
The lattice filter is another optimal least-squares 
predictor which can be applied to the linear deconvolution 
problem. The lattice, or ladder filter derives its name 
from the cascade form of its signal flow graph. Basically, 
the lattice filter provides an alternative to the 
transversal filter for modeling a signal. It can be viewed 
as a Gram-Schmidt orthogonalization of the incoming data. 
This will be elaborated upon in subsequent paragraphs. To 
date, much of the work done with lattice filters has been in 
the areas of speech analysis/synthesis, seismology, and in 
high-resolution spectral estimation [{Ref. 7:p. 841]. Prior 
to developing the lattice filter equations, key mathematical 
concepts which apply to this discussion will be reviewed. 
2. Mathematical Background 
A central concept behind the development of the 


lattice filter is that of orthogonality. Two vectors are 
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orthogonal if their inner product is zero. The geometric 
interpretation of this is that the vectors are at right 
angles to one another. As an example, if the random varia- 
bles uandv are the basis vectors of a two-dimensional 
space, they are orthogonal if and only if their inner 
product <u,v> = Efuv] = 0. In the case where the linear 
Space is spanned by the random variables of a data vector, 
two vectors are orthogonal if the corresponding random 
variables are uncorrelated and if one or both have zero 
mean [Ref. 8:p. 92]. 

A related theorem is the orthogonal decomposition 
theorem, which states that any random variable may be 
decomposed uniquely with respect to a subspace S into two 
mutually orthogonal parts, one part which is parallel to §S 
(i.e., lies in S) and the other part orthogonal to S. That 
is, y can be written y = v + e where e is orthogonal to y 
and to the basis vectors which span the subspace 5S, and 


where o is defined by a linear combination of the basis 


vectors. lf S is spanned by the basis vectors 
{u(1),u(2),...,u(M)}, then ¥ = y a(i)u(i) where it can be 
shown that the coefficients ae given by the equation 
eae) = Te art a trefo lep. 13] 


The orthogonal projection theorem adds to this’ by 
Stating: the orthogonal projection of a vector y onto a 


linear subspace S is that vector in S which lies closest to 


39 


y with respect to the distance induced by the inner product 
of the vector space [Ref. l:p. 15]. Restated, this simply 
means that 7 is the best estimate of y that can be made by a 
linear combination of the basis vectors of S in a minimum 
mean-squared error sense. Figure 2.5 shows the projection 
of y into the subspace S. 

Another important concept in developing the lattice 
filter is that of Gram-Schmidt orthogonalization. The Gram- 
Schmidt procedure is a recursive orthogonalization process 


which generates a set of mutually orthogonal basis vectors 


(all eC Ze from a given set of basis vectors 
fey (ee eZ) eee eve Millon The procedure is initialized by 
Tetting ull) = JyCL). Next, y(2) is decomposed with respect 
TOne ul das. That part of y(2) which is orthogonal to u(1) 
becomes u(2). Next, y(3) is decomposed with respect to the 


subspace spanned by {u(1),u(2)}. That part of y(3) which is 
orthogonal to this subspace becomes u(3). In general, the 


new set of orthogonal basis vectors is defined by 


n-l 
u(n) = y(n) - >) b(n,i)u(i) (2.56) 
ite 
2 -1 
for 2<« n < M and where b(n,i) = Ely(nju( py bluawe 
Using the orthogonal decomposition theoren, thas is 
equivalent to y(n) = y(n) + u(n) where y(n) is the best 


estimate of the n-th component of the data vector based on 


the previous n-l components and u(n) represents the 
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Figure 2.5 Projection of y onto Space S 
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he 


estimation error. If b(n,n) = 1 iS Introducedy eur. 
3 
y(n) = dy, b(n, i)u(i) for 1 < sna te (2. 5 dam 
eo 


This can be written in matrix form as 


y = Bu- where (2.5 
T 
y = Ly(ljigmieGdle - ay , 
T 
u-= [ut l),9ut 205 2 eee ’ 
é 
1 0 0 eo ¢ 66@ 0 
bi'Z;. lt) 1 0 : 0 
a — bt(3, 1) bias 2) 1 os 0 
bie; £) b( M52} i 
This is a convenient notation for representing the 


transformation from a set of correlated basis vectors y to 
an uncorrelated set of basis vectors u. The bases y and u 
Span the same M-dimensional subspace S, but there are no 
redundant correlations between the basis vectors of set us. 
Since the basis of u is uncorrelated, its components) u(i1) 
(1 = 1,2,...,M) are referred to as innovations because’ each 
additional component contributes completely new information 


to the estimate of y. [Ref. l:pp. 16-18] 


Finally, the transformation y = Bu corresponds to a 
LU (lower upper)-Cholesky factorization of the correlation 
i 
Matyix ot wy. By definition R = Elyy J]. Substituting wien 
yo 
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y results in 


Tae T 
R = BE{uu JB = BRB O25. 58 } 
yy uu 
T 
where B is lower triangular, B upper triangular, and R is 
uu 
a diagonal matrix since the basis vectors) u(i) are 
uncorrelated with one another. 
3. Derivation of Lattice Filter Equations 
Now the lattice filter order update equations will 


be developed. A random signal 


output of a causal, stable, 
driven by a stationary, uncorrelated, 
{e(l1),e(2),...,e(P)} [Ref. eo 


autoregressive model is defined by H(z) = 


-1 -2 
A(z) = +a(Z)z 


P 


lta(l1l)z 


The filter A (z) has several names: 
Pp 


inverse filter, or analysis filter. 


written as 


Pp 
ee c(h) = >, a(i)y(n-i) 
i=1 


If the predictor y(n) 1s introduced, 


y(n) - y(n) = e(n) 


Pp 
- ). alidy(n-i). 


=] 


where y(n) = Now, 


y(n) at time (n-1) based on the previous P samples, 
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y(n) can be modeled as 


linear filter H(z) 


+...ta(P)z 


the 


which is 


white noise sequence 
309% A P-th order 
1/A (z) where 
Pp 
-P 


(ozo oO) 


prediction error filter, 


The signal y(n) can be 
(2 760:) 

this becomes 
(2 Sorls) 


7(n) is the estimate of 


iy, 


Ih 
ip 
cy 


vin-P+1),..., yvin-l1)} and e(n) is equivalent to the estima- 
tion error. The estimate Fin) can be thought of as’ the 
projection of the (P+1)-dimensional y(n) vector onto the P- 
dimensional subspace spanned by the components of the data 
vector Y = leet). sin In general, the 
past values of y(n) are correlated with one another. There- 
fore, the random variable components of Y do not generally 
form a set of mutually orthogonal basis vectors. The Gran- 
Schmidt orthogonalization procedure removes these 
correlations and transforms Y into a eset of mutually 
orthogonal basis vectors which span the same subspace. 
Additionally, the predictor coefficients a(i) are replaced 
by the reflection coefficients K . This results in the 


1 
lattice filter. 


In order to obtain the optimal estimator, the 
predictor coefficients which minimize the mean-squared 
prediction error eis cue must be determined. This problem 
was discussed in the introduction to this chapter. The 


resulting matrix form of the normal equations is 


RCO) RUD) GRGZ) een eae) if Ep 
Rei} RO) Ro) ere. a(1) 0 
R62 RCA) RCO r ee canes a(2) = 0 (2.62e 
R(P) R(Q) a(P) 0 
where R(k) = Efy(ntk)y(n)]}] and E is the minimized mean- 
P 
squared prediction error for the P-th order filter. The 
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sample autocorrelation matrix components are calculated from 


the length N data sequence by 
N-1-k 
R(k) = (1/N) )  y(nt+k)y(n) f2eeoD 
n=0 


for 0 < k < P. and where P < N-1l. [Ref. l:p. 150] 


Now there are P+l equations and P+l unknown model 


parameters {a(1l),a(2),...,a(P);E }. The P+l equations can be 
P 

solved by inverting the autocorrelation matrix. This 
3 Z 

requires O(P ) operations and O(P ) storage locations. The 


P+1 equations can be solved more efficiently by taking 
advantage of the matrix’s Toeplitz structure by using 
Levinson’s algorithm. Levinson’s algorithm reduces’ the 
required number of operations and storage locations to oa. 
and O(P), respectively [Ref. l:pp. 150-151]. Being a 
recursive procedure, Levinson’s algorithm permits’ the 
calculation of the (P+l)-st order model parameters by using 


the previously determined P-th order model parameters. The 


matrix form of the algorithm is given by 


45 


(i 
he 
n 


a (1) a (le) a (P) 
P+1 iE P 
a (2) a (2) a (P-1) 
P+1 P P 
: = : amen : (2.049 
‘ : P+1 ; 
a (P) a (P) a (1) 
P+1l P P 
a (P+1) | 0 | 1 
P+1 | 
where 
P 
). a (i)R(P+1-i) 
V=0> =P 
K = (2.65m 
P+1 P 
». a (i)R(i) 
LS0 P 


st 1 


The P subscript on the "a" parameters specifies the P-th 


order prediction error filter A (zZ2) Whereas the index 
identifies the appropriate term i the A (Zz) polynomial. 
K is called the (P+1l)-st order eee or PARCOR 
ioe cia correlation) coefficient. The PARCOR coefficient 


K represents the true, or direct, correlation between 
ei and y(n) with the effects of the intermediate 
variables (i.e., y(n-P),y(n-Ptl),...,y(n-1)) removed. The 
recursive form of Levinson’s algorithm is written as 
a (m)= a (m) —- K a (P+l-m) for 1<m<P, (2.66a) 
P+1 iE Pa ie 


and 
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a (P+1) = -K for m=P+l. 20.6.) 
P+1 P+1l 


This shows that the reflection coefficient for each stage P 
is equal to the highest coefficient of A (za). There is a 
one-to-one correspondence between the ak coefficients 
and the coefficients of the transfer function A (z). The 
transfer function or, equivalently, the a ae 


matrix R = Efy(n)y (n)] uniquely determine the reflection 


coefficients [{Ref. 7:p. 829]. Taking the z-transform of 


equation (2.66) results in 
A (z) = A (Zz) - K B (z) (-2ea6 7 2 ) 
P+1 P P+1l1 P 
where 
-1 -2 -P 
A (z) =lta (1)2 ta (2)zZ +...t+a (P)z ; (2.675) 
P P P P 
and 
-P -1 
Beez ee nz A 2), (2.67c) 
P P 
-1 -2 -(P-1) -P 
= a (P)ta (P-1)zZ ta (P-2)z2 +...+a (1)z +Z 
P P P P 


Due to the effective folding about the axis and the shifting 


to the right of the sequence A (z), the polynomial B (za) is 
P P 
called the reverse of A (2). Taking the reverse of both 
P 
sides of equation (2.67) yields 
-1 
B CZoee= 2 Bb CZ )-=— K A (2) (2.606) 
P+1 P P+1 P 
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Combining equations (2.67) and (2.68) into matrix forme ie 


-1 
A Cz) 1 -K Z A (2) 
P+1 P+1 P 
= -] (2.69) 
B (z) -K Z B (z) 
Peta P+i1 Pp : 

The forward prediction error associated with 
predicting y(n) from the previous P samples {y(n-P),...y(n- 
1)} is written as e (n) = y(n) - y (n) where the subscript P 

P PE 
represents the filter order. In z-domain notation, this 
becomes 

Be (2) A ez) (2.7m 

P P 
Now, the backward prediction error r (n) is introduced. It 
P 
is defined as the error in predicting ( or actually 


smoothing) y(n-P) from the future P- samples {y(n-P+l,..., 


vin) s It 1s written as 
r(n) = y(n-P)ta (1l)y(n-P+1)+...ta (P)y(n) (2.¢la) 
Pp P P 


or in z-domain notation as 


E (z) = B (z)Y¥(z) . (2. 7iuep 
P P 
The optimal backward predictor coefficients minimize the 
Z 
mean-squared smoothing error E[r (n)]. Since the optimal 
P 


backward and forward predictor coefficients are mirror 
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2 Z 
images of each other, unineme ir (n)! = Efe (n)|J.. That is, 


the backward and forward prediction error vectors have the 
same norms [{Ref. 8:p. 101]. 

As will be shown, at each instant n, the backward 
prediction errors are mutually orthogonal (i.e., uncorre- 
lated if assumed to be zero mean) {fRef. ps ie Oe. 
Therefore, Tete 1s the backward prediction errors, 
mene ,p-0,1,...,M, where M is the filter order, which form 
a new set of basis vectors. To demonstrate that the 
backward prediction errors are mutually orthogonal, let M=3 


and write the corresponding equations for P = 0,1,2,3 in 


matrix form: 


ONG 


r (n) 1 0 0 O y(n) 
0 
meetin) a (1) li 0 0 y(n-1) 
1 1 
= G2. 72a) 
r (n) Ai 2 joe | Le) 1 0 y(n-2) 
Z 2 Zz 
mT) aumesor a (2) “ant hyeo y(n=3) 
3 3 3 3 ‘ 
Gi) ea iy ( 1) : CZ Zp) 


Note that the first column of L contains the negatives) of 
all the reflection coefficients. Now examine the covariance 


matrix of r(n): 
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T i: 
R = E[r(n)r (n)] = E(LyY( mee byte ( 22720 
tots 
T i T 
= LElyi(nj)y (n) Jb "= Eke oe 
yy 


Since the normal equations are satisfied within the matrix 


products above, R reduces to a diagonal matrix which 
tone 
verifies that the components of r(n) are uncorrelated. That 


is 


i 
R = LR L = D = diag{i eee (2.4 
rr Vy O- SI ees 


where E is the minimum value of the mean-squared prediction 
P 2 
error which is given by Ef[e (na) “s]'% It can be shown that 
2 P+1 
E = (1-K JE where the recursion is initialized with 
P+1l P+1 P 2 
the value given by E = R (0) = Ef[y (n)] {Ref. 8: p. 10pm 
0 yy 
Therefore, the prediction error decreases by a factor of 


Z 
(1-K ) from one lattice stage to the next. Now, since the 


remit of r(n) are uncorrelated, equation (2.72) 1s 
equivalent to the Gram-Schmidt orthogonalization of the data 
vector y(n). Also note that rewriting equation (2.74) as 
R = , ae corresponds to a LU-Cholesky factorization of 
She covariance matrix of y(n). 

To complete the derivation of the lattice recursion 


equations, multiply both sides of equation (2.69) by Y(z) to 


get 
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E (z) 1 -K z, Eire 


P+1 P+1 P 
om = -1 my 2 1D) 
E (z) -K Z, EK (2) 

P+1 P+l P 


These equations are transformed into -+tthe time domain to 


obtain the final result: 


e (n) e (n) - K r .(n-1) (2 Oa) 
P+1 P P+1 P 


rar) =r tn-1) - Ke (n) . (2.76b) 
P+1 P P+1 P 


These equations, which are recursive in both time and order, 
define the signal flow graph of the analysis lattice filter, 
shown in Figure 2.6. For a given time instant n, the equa- 
tions are eevee etn sively in order. The inputs into 
the first stage of the lattice filter are e ({n) = r ({n) = 
cat) . The successive stages of the er etetes the 
Successively higher order forward and backward prediction 
errors. The output from the final stage yields the desired 
M-th order forward prediction error, while all lower order 
prediction errors are available at intermediate stages. 
Since the backward errors are generated from the y data 
vector, the analysis lattice filter actually implements 
the Gram-Schmidt orthogonalization; all that is required to 


implement the lattice are the reflection factors 


{K »K emery} 
i 2 M 


Ol 





Figure 2.6 Analysis Lattice Filter 


2 


Equation (2.76) can be manipulated to obtain the 


it tice form of the synthesis filter H(z) = 1/A (2): 
M 
e (n) =e (n) + K r (n-1l) iZel tad 
P P+1 P+1 P 
i (n) = r (n-1) - K e (n) . (2 Cb) 
P+1 iE P+l1 P 
Figure 2.7 shows the corresponding signal flow graph. When 


the input to the synthesis filter is the forward prediction 
error sequence e (n), the output is the original sequence 
mam) . in nie oh the synthesis filter to be stable and 
causal, all M zeros of the prediction-error filter A (z) 
must lie within the unit circle in the complex em. A 
necessary and sufficient condition for all the zeros to be 
inside the unit circle is that the magnitude Be enh of the 
reflection coefficients {K ,K ,..-,K } be less than one. 
ec M 

[Ref. l:ipp. 168-169] 

The reflection factors can be evaluated by several 
methods. The various methods arise due to different 
definitions of the optimality criterion. The criterion 


used here of minimizing the mean-Squared prediction errors 


leads to what MAKHOUL calls the forward and backward methods 


[Ref. 9]. They are, respectively: 
2 
K =eihre in) (n=l) ) 7 Eltr (n=) ) (2.78a) 
P+l,e P P Pp 
2 
K = hre (nj)r (n—-l)]) / Efe ({(n) } z (2Ze1Sb) 
P+l,r P P P 
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Figure 2.7 Synthesis Lattice Filter 
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a 2 


Assuming stationarity, and since E[r (n) ] = Efe (n) ] as 
P P 
previously discussed, then the forward and backward 


reflection coefficients are equal. That is K = = 
P+1 P+l,e 


K »- The Schwarz inequality implies that ;K ; < 1 for 
eri ,r P 
each P=1,2,...,M [Ref. 8:p. 104]. An alternate technique for 
calculating the reflection coefficients is the geometric 
mean method. It was introduced by ITAKURA and SAITO in 
their work of developing a digital filter structure for 
time-domain speech analysis [{Ref. 10]. The corresponding 
PARCOR coefficient is given by 
Efe (n)r (n-1)] 


Pr P 
K Se CZ oO ) 


P+1 2 2 1/2 
Penryn — 1) 
P P 
These PARCOR coefficients are guarenteed to have magnitudes 
less than one [{Ref. 7:p. 840]. 

A third method was used by BURG in the maximum- 
entropy method of spectral estimation [Ref. 11]. This is 
the harmonic-mean method. The harmonic-mean method seeks’ to 
minimize the sum of the forward and backward prediction 
error variances, Efe + E[r eee Minimizing this 


P P 
Sum with respect to the reflection coefficients results in 
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Mi 


FREFROUUCED Al GUVEKINMENt BEAreciIVNoOc 


2E[e (n)r (n-1)] 


P P 

K ee (2.008 

P+1 Za 2 

Efe (n) + E{r (n-1)] 

iB p 

Again, the Schwarz inequality verifies that K has magnitude 
P 
less than one, guaranteeing that the synthesis filter is 
causal and stable [Ref. l:p. 189]. 
a The Generalized (Analysis) Lattice Filter 


The preceding development of the analysis’ lattice 
filter equations assumed that the data sequence {y(n)} 
represented a time sequence with stationary statistics. In 
this section, a more general linear prediction problem 
will be considered. No special assumptions) are made 
concerning the data. The data values need not be delayed 
versions of each other; they need not even represent a time 
sequence. This is the approach taken by LENK in developing 
the generalized order update equations [Ref. 12:p. 85]. The 
resulting generalized form of the lattice filter makes it 
suitable for multidimensional and nonlinear Signal 
processing applications. 

Definitions associated with the normalized form of 


the generalized lattice filter will now be introduced. 


First, the components of the length (M+1) column vector 
ty? ie representing a single realization of the random 
process Y, are designated by ya : A= O,1,...,M. The forwawd 


prediction error in estimating the (ntl)-st element from 


06 


the previous N elements of the data vector is written as 


N M N x 
2 = y h (ntl)y (20 1) 
n+l x= 0 A 
where the length (Mt1) row vector of coefficients is 
given bv 
N N N N 
{h Gite 1=(.0),.5.,0,—h ,-h 9 ah leone Ol a ( 2<Oe ) 
d n-N+l n-N+2Z n 


The backwards prediction error associated with predicting 


n-N 

Vv from the next N elements of the data vector is given by 
N Mae N x 
r = ) h,(n-N)y (283) 
n-N A=0 x 


where the associated coefficients are given by the vector 


N N N 


~ ~ rw 


N 
{h (n—N) 1=(0,...0,1,=h ,-h iets Opes sO le (2.64) 
Xr isle ar I oiler n 


The norm of the forward prediction error is defined as 
N Ni pce 2 
lle Ud = (E{(e ) }] ° (2285) 
n+1 n+1 
The norm of the backward prediction error is defined in a 


Similar manner. Now, the normalized forward and backward 


N-th order prediction errors are defined by 


N N N M N x 
e =e {/ ‘tle ices 3 a (ntl)y ; (2.86a) 
n+1 n+] n+l A790 » 
and 
_N N N M N x 
r =i (fea (Ts oa {j = bl.(n-N)y (2.86b) 
n-N n-N n-N =0 A 


om | : 


where the normalized prediction error coefficients are 
defined as 
N N N 
a (ntl) =] ne (ntl) eile eae (2.87a) 
r A her 
and 
f 
; N _N N 
b.(n-N) = h.(n-N) / flr | : (2.87b) 
’ x » n-N 
[Ref. 12:pp. 85-87] 
Using the normalized form of the generalized 
Levinson algorithm, LENK demonstrated that equations (2.87) 
could be updated recursively through the relation 
i$ 
N+1 N 
i a. (n+l) n a. (nt1) 
ni ‘ = Gi lee (2.88) 
ae N+1 N+1 N 
Nios = -N 
cue Bey a ah) b(n-N) 
where the partial correlation (PARCOR) coefficient is 
n N N 
K = eee r } ; (2.89) 
N+1 n+l n-N 
and where 
n 
n 1 il ~ K 
9 (xk een N+1 (2.90) 
N+1 n 2 n 
1 - (K ) - K 1 
N+1 N+1 ; 


The recursion is started for order N=0 with the initial pre- 


diction error values fore 0, eee 


8 


~-»M given by the vectors 


ntl 
oie) yet Oe ery Oy I I y bi} ,0,...,0] (2779 Wan) 
n 


Poe (moles [0,...,0, 1/7 | ty Wee cueny Onl eee (Zao 1b) 


+» oO + © 


Multiplying both sides of equation (2.90) by the data vector 


ty? | yields the desired error order update equations: 


N+1 N 
e e 
nt+1l n n+l (2.92) 
- §(K ) 
N+1 N+1 N 
r rT 
| n-N | | n-N | 


The corresponding signal flow graph for a single lattice 
filter section is shown in Figure 2.8. Figure 2.9 depicts a 
third order generalized lattice filter. [{Ref. 12:pp. 94-99] 
LENK then proved that the reflection coefficients 
(i.e., PARCOR coefficients) could be calculated directly 
from the data vector’s correlation matrix by utilizing the 
generalized Schur algorithm. In LENK’s notation, the 


A 
RP = 


components of the correlation matrix are written as 


E{yP oo }. Then the reflection coefficients are given by 
n-N 
OY Grit) 
n N 
K ee : (229:3)) 
N+1 n-N 
GL (n-N) 
N 
where 
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b. Alternate Representation 


Figure 2.8 Generalized Lattice Filter 
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Figure 2.9 3-rd Order Generalized Lattice Filter 
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MEPrIMUUULCED AL GUVEMNiCint EAtoie 


- 


N A 


a M Y 
Q (ntl) = E(@ y} = fa (mtl)R , (2.94a) 
N+1 n+l fiat 
and 
x Nae M N pA 
GL (n-N) = E(r  y } =e ib euoeny hae (2.94b) 
N+1 n-N } =0 


Equations (2.94) can be updated through the recursion 


relation 

» IN 

(ntl) CO (ntl) 
N+1 n N 

= (K ) (2.955 

A N+] » 

L (n-N) iG (n-N) 
N+1 | N 


To start the recursion at order N=0, the values of the 


parameters for A=0,1,...,M are given by the vectors 
Xx n+l -1l (n+1)O0 (nt+1)1 (n+1)M 
(Oe ne ya ty, {| [R »R ees: ] (2.96a) 
0 
DN hey) ok nQ nl nM 
(OC myl=iiy ti 1k OR 7 ee (2.96b) 
0 


[Ref. 12:pp. 100-101] 
9. Lattice Filter Advantages 
The lattice filter has several advantages compared 
to the direct, oor transversal, form of the prediction error 
filter Aine First “of, sal due to the built-in 
M 


orthogonalization incorporated into the lattice, successive 


lattice stages are decoupled. Therefore, the reflection 
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Peetticients K , P=1,2,...,M are independent of the filter’s 
minal order. sar M-th order least squares prediction filter 
contains all the prediction error filters of lower order. 
Restated, the first P sections of the M-th order filter form 
the P-th order prediction filter. Therefore, lattice stages 
may be added or subtracted from the existing lattice filter 
without having to recalculate the already determined 
reflection coefficients. In contrast, when the order of the 
transversal filter 1s changed, allthe A (z) prediction 
error filter coefficients must _ be ee ct ca. AS a 
consequence, to specify all filter orders up to M requires 
storage of M(M+1)/2 predictor coefficients, while only M 
reflection coefficients need to be stored. (Ref. 7:p. 830] 

Another advantage of the lattice over the 
transversal filter is that the output of each lattice stage 
is a least-squares prediction error. Therefore, if the 
desired order of the predictor is not known in advance, the 
output of the various stages can be monitored to determine 
what filter order is adequate. [Ref. 8:p. 106] 

Due to the quantization inherent in implementing 
Gigital filters, round off noise is introduced.: Studies 
have shown that the performance of the lattice filter is 
far superior to that of transversl filters for finite word 


length computations. Furthermore, the filter zeros are less 


sensitive to quantization errors in the reflection 
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coefficients than for quantization errors in the transversal 
filter coefficients. Finally, Since the magnitude of the 
reflection coefficients is less than one, it simplifies the 
task of establishing the quantizer overload point. [Ref. 
Sem ee 

Another advantage the lattice filter holds over the 
transveral filter is that the lattice filter is minimum 
phase (i.e., all its zeros are inside the unit circle _ so 
that the synthesis model is stable) if and only if the 
magnitude of the reflection coefficients is less than one. 
There is no comparable test for the transversal filter 
coefficients to determine whether the synthesis model is 
stable. [Ref. 8:p. 106] 

6. Linear Lattice Filter Applied to Deconvolution 

The transfer function of the lattice filter is 
determined by the reflection coefficients. In turn, the 
values of these reflection coefficients are uniquely 
determined by the prediction error filter transfer function 
A (z), or equivalently by the autocorrelation sequence R(k), 
oo O,1,....M (Ref. T7:p. 829]. Similarly, equations (2. 
and (2.71b) indicate that the transfer -functions from the 


input y(n) to the outputs e(n) and r (n) are A (z) and 


M M M 
B (z), respectively. Therefore, the analysis lattice filter 
M 
is equivalent to the whitening filter A (z); that is, it is 
M 
the inverse of the system model H(z) = 1/A (z). This is the 
M 
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basis for using the lattice filter ina deconvolution 
application. 

In order to recover the input signal x(n) from. the 
system’s output signal y(n), it is necessary to determine 
the inverse of the system transfer function H(z). The 
initial step is to conduct an "identification experiment" to 
estimate the system’s parameters (either the autoregressive 
filter coefficients or the latence filter reflection 
coefficients). To accurately estimate the parameters, the 
chosen experimental input signal must be sufficiently rich 
in frequency content, or more formally, persistently 
exciting so that it excites all the modes of the system. A 
sequence is said to be persistently exciting of order n if 
its (nxn) autocorrelation matrix is nonsingular ([Ref. 
13:pp. 70-71]. By using the techniques presented in section 
D.3 of this chapter, the resulting output sequence y(n) can 
be processed to yield the desired reflection coefficients. 
When the analysis lattice filter is implemented with these 
coefficients embedded in its structure, the lattice becomes 
the inverse of A (z) = 1/H(z). This procedure is depicted 
MreFigure 2.10. . 

The application of linear lattice filters to decon- 
volution was simulated using the FORTRAN programs’ LININV, 
ATOCOR, LEVIN, AND LATICE. These programs are provided in 


Appendix A. The simulation was conducted for a_ second 
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Figure 2.10 Deconvolution Using Lattice Filter 
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order, autoregressive, LTI system defined by the equation 


y(n) = x(n) - (0.6y(n-1) + 0.08y(n-2)) . As previously 
mentioned, modeling the system was the first step in the 
deconvolution process. In order to identify the "unknown" 
parameters, the system was excited by a zero mean, unity 
variance, Gaussian white noise sequence. The output 
sequence y(n) was processed by subroutine ATOCOR ZO 


calculate the components of the autocorrelation matrix Re. 
Then subroutine LEVIN implemented Levinson’s algorithm - 
evaluate both the prediction error filter coefficients 
and the associated reflection coefficients from the autocor- 
relation matrix. The actual output of subroutine LEVIN is 
the lower triangular L matrix described in equation (2.72); 
the i-th row of L contains the i-th order prediction error 
filter coefficients listed in reverse order, and the first 
column contains the negative of the reflection coefficients. 
Once the reflection coefficients corresponding to 1/H(z) 
were calculated, they were embedded in subroutine LATICE 
which implemented the analysis lattice filter equations 
(2.16). With the inverse filter of H(z) now available 
(1.e., the analysis lattice filter), H(z) was driven by an 
"unknown" signal x(n). The output of H(z) was then fed into 
the lattice filter. An approximation to x(n) was recovered 


at the forward error output signal from the last stage of 


the lattice. 
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Simulations were run both with and without’ the 
presence of measurement noise. Table 2.1 displays’ the 
resulting L matrices for one simulation. The case of no 
measurement noise 1S shown in Figures 2.11 through 2.14. 
Figure 2.11 is a plot of the input signal x(n), Figure 2.12 
depicts the system output y(n), and Figure 2.13 shows’ the 
lattice filter output ee As can be seen, the results of 
the inverse filtering were excellent--the x and & curves are 
identical. This is verified by Figure 2.14 which shows that 
the mean-square error between x(n) and x(n) 1s nearly zero. 
The running average mean-square error was calculated using 


the equation 





2 


n 
MSE(n) = (1/n) ¥ (x(a) - x(i)) (2.97) 
= 1 





The simulation was then repeated for a nonzero 
measurement noise. Here, the added measurement noise was a 
zero mean, 0.0025 variance, Gaussian white noise sequence. 
Using the same input as for the previous case, the outputs 
of the system and lattice filter are shown in Figures 2.15 
and 2.16, respectively. Figure 2.17 is a plot of the mean- 
Square error between x(n) and &(n). The results in Table 2.1 


show that even with the presence of measurement noise, the 


system parameters were still accurately identified. The 
lattice filter recovered the basic shape of x(n), however, 
Ue could not remove the additive measurement noise. 
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Therefore, the output of the lattice filter is a noisy 


or "fuzzy" version of x(n). Additional filtering is required 
to remove the noise component of x(n). The results 
presented here are representative of those obtained for 
Simulations involving other’ stable, autoregressive, LTI 
systems. 
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TABLE 2.1 


RESULTS OF MODELING A LINEAR SYSTEM 


Modeling Problem: 


-1 -2 


System: H(z) = 1/{1 + 0.62 + .08zZ J 


System input: Gaussian white noise, N(0,1) 


Number of white noise realizations used: Zo 


a. 


Number of points per realization: 5; 060 


Modeling Results: 


No Measurement Noise 


1 0 0 
ie O05 LO 0 
0.07% O2598 1 
Order ,P Ep Kp 
0 1.445 -- 
il Oe dg0 -0.555 
Z O23 93 -0.077 


Measurement Noise is N(0,0.0025) 


1 0 0 
lL = U2 000 Ie 0 
OSG <O<59S 1 
Order ,P Ep Kp 
0 1.449 -- 
1 1.007 =0 5559 
Z 0.996 -0.077 
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Figure 2.11 System Input, x(n) 
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Figure 2.12 System Output, y(n) 
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Figure 2.13 Lattice Filter Output, x(n) 
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Figure 2.15 System Output Plus Measurement Noise 
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Figure 2.16 Lattice Filter Output, x(n) 
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III. NONLINEAR DECONVOLUTION 


A. INTRODUCTION TO NONLINEAR SYSTEMS 

While the previous chapter dealt solely with linear 
systems or plants, this chapter will address the inverse 
filtering problem involving nonlinear systems. The chapter 
starts with a brief introduction to modeling nonlinear 
systems. Then it quickly proceeds to extend the generalized 
linear lattice filter results of section II.D.4 to a 
nonlinear analysis lattice filter. The nonlinear lattice 
filter is discussed in detail. To conclude, results from 
numerous’) simulations involving the lattice in deconvolution 
applications are presented. 

System linearity is defined in terms of the principle of 
superposition. If the rule by which the system transforms 
the input x(k) into the output y(k) is represented by the 


operator T, then the system is said to be linear if and only 


Tete 
TlLax (k) + bx (k)] = aT[x (k)] + bT[x (k)] (3a 
1 2 1 2 
= ay (k) + by (k) 
1 2 
for arbitrary constants a and b. The system is7~ nonlinear 
if equation (3.1) is not satisfied. Estimation of the 


parameters of a nonlinear system is a complex problem. 
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One approach to the parameter estimation problem for 


nonlinear filters is the Bayesian approach. This approach 


leads to various approximate solutions for the parameter 
estimates. In the Bayesian approach, the parameters’ are 
considered to be random variables. The parameter vector h 
is considered a random vector with a probability density 
function p(h). Introducing the observation vector y(t) and 
the input vector x(t), both of which contain data up to 
time t, the a posteriori probability density function for h 
mempihty(t),x(t)). One possible choice for the estimate of 
the h vector is to select the conditional mean nee = 
Efhly(t),x(t)J]. Selecting this as the estimate minimizes the 
variance of the parameter estimation error. Another 
possible choice is to select the nen which maximizes 
mil y(t),x(t)). This most likely value is known as’ the 
maximum a posteriori (MAP) estimate. Finally, the Bayesian 
approach also leads to the extended Kalman filter for 
nonlinear state estimation problems. (Ref. 13:pp. 32-41] 
Another popular method for modeling nonlinear systems is 
based on the Volterra series. The series is named after the 


mathematician Vito VOLTERRA. The first person to apply the 


series to nonlinear systems was Norbert WIENER. If the 
"black box" approach is taken towards a nonlinear, time- 
invariant system, the relationship between the input x(t) 


and the output y(t) can be represented by the Volterra 


ao 


ee FO FY LF 


ee ee er ew - em ow 
aaa Se Se Pe FY Y, 


-™ 


zc? 


series: 


oO 
y(t) - fr (T ) xu a? 
aoe 1 i 
co 
lf h (T ,T )x(t-T )x(t=7T aie 
wea 2 1 2 nae 
CO 
f h (T ,T ,T )x(t-T )x(t-T )x(t-T )dT dieu 


+ 
WE tte 1 2 3 Te 


Te 6 ¢ 6 


L ft, h —- pooosgl )X(t=—T ex (toe eee 
n 1 n 1 n 


(3 92a 


where n =1,2,... and -h (T ,...,1T ) = 0 £f6r any T =e 
n | n J 
j = 1,2,...,n . The functions h (T ,...,1f ) are CGaiiuee 
no n 
Volterra kernels of the system. This equation is a 
functional series. That is, it performs an operation on the 
function x(t) which results ina number for y(t). If the 


n-th order Volterra operator, H , is introduced where 
n 


y (t) = H [x(t)] (3.3) 


n 
co CO 
:| a h (T ,.-+,T )}xX(t—-T ) 3.x ( t—T See eee 
n 1 n Il n il n 


then the series can be expressed in operator notation as 


y(t) = H [x(t)] + H [x(t)] + .22) + HH o[x(t) ] eee (3.4) 
1 Ze n 
co co 
= ). H (x(t)] =") yee 
he) aon n= gn 
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Figure 3.1 is a graphic representaion of this equation. 
[Ref. 14:pp. 7-9] The H operator is a linear operator, 
H is a quadratic estate a and, in general, the H operator 
es ves the term x(t) to the n-th power. i‘ 
B. GENERALIZED NONLINEAR LATTICE FILTER 
1. Introduction 

In this section, ie will be demonstrated how the 
generalized lattice filter introduced in section I1I.D.4 can 
be applied to modeling nonlinear systems. The development 
will start by looking at the discrete form of the Volterra 
series. Based upon this series an alternate tensor notation 
representation will be introduced. The results of section 
II.D.4 will then be extended to handle a_ two-dimensional 
field of data which will result in the generalized nonlinear 
lattice filter. 

2. Nonlinear Lattice Filter Development 

The discrete form of the Volterra series of equation 

(3.2) is given by 
y(k) = h + h (n )x(k-n ) C355) 
0 1 1] j 
+ h (n ,n )x(k-n )x(k-n ) +... 
Zuele e2 1 2 


Using LENK’s tensor notation, an equivalent form of equation 


(3.5) is given by 
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Volterra Series Representation 


Figure 3.1 
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A\ A Az 
eee + 


Vues xX ~~ H .xX xX » where (3.6a) 
A A\Az, 
ri a 
x (k) = [x(k),x(k-1),...,x(K-Ar),-..,X(K-N) ] (3.00) 
meen = 6 0, 1,4. wy N. As an example, if N=l and equation 
(3.6a) is truncated to the three terms shown above, then 
this equation can be rewritten as 
y(k) = H + H x(k) + H x(k-1) + H x(k)x(k) 
0 i: 00 
+H x(k)x(k-1) + H x(k-1)x(k) 
O1 10 
+ H x(kK-1)x(k-1) . (3.7) 
11 
Based on this tensor notation, LENK introduced an alternate 
representation for the nonlinear’ system. Instead of 


defining the components of the vector x (k) as the present 
and past values of the signal x(k) as in equation (3.6b), 


the components are redefined in terms of x(k) raised to the 


A, power ford; =O) ling Atrio Neo 0k NieeekssS 
r; 0 1 2 N ih 
x (ie = (x Ok), x Blin, x Ce. x Ck) ] ( Sree) 
2 N T 
=e ON eee 1) Gee KO ) | 
Lor nr; = 0,1,2,...,N. Now the output of the nonlinear system 


is defined by 


y(k) = ae aac N)H (3.9) 
ee ae : Ne ; 
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for Dee = 0,1,2,...,;P. The variable P 2ivecs@upee. 
of the filter and N defines the filter’s finite memory. The 
term ay plays a role similar to that of the Volterra 

o N 
kernel; it can be considered a (P+1l)-order tensor. To 
clarifly the meaning of this notation, the following example 
with N=1 and P=2 is provided: 

i) = Hee eee 
y(k) = A> ° 
- H +H x(kK)+H x(k-1)4+H =x(k)x(k-1)+H x (k) 
00 10 O1 i] 20 
2 Z Z 
+ H x (k-1)+H x (kK)x(k-1)+H x(k)x (k-1) 
02 Dal 2 
2 Z 
oe 
[Ref. 12:pp. 39-48] 

The above nonlinear model for y(k) iS a moving 
average (MA) model: the output is defined in terms of the 
input Signal. The next step is to establish an 
autoregressive (AR) model where y(k) is estimated in terms 
of its past values. This type of model is useful when the 


input signal is not readily available. An autoregressive 


model is obtained by redefining the observation vector in 


terms of y(k) vice x(k). Then the AR model is given by 
A\ AN 
y(k) = y (k-1).. 3k ND eee (3. ie 
A, An 
fie t A; = eZ ga ages If the system is driven by a_e white 
noise signal u(k), and if the system’s parameters) are 
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An 
approximated by H then the estimation error is 


ee 


At). . ey NEN) H ee i 


e(k)=y(k)-y(k)=Ly aD 
\ N 


A AN iN 
- Cy Us-1)..+y (K-N)H) + .-Q (3.12) 


If the system parameters are known exactly, then the estima- 
tion error and the white noise sequences are equal--that is 
e(k)=u(k). One note concerning the nonlinear AR model: . 
unlike the linear AR model whose stability is easily 
determined (i.e., the system is stable if all its poles lie 
within the unit circle in the complex z-plane), it is 
difficult to judge for which class of inputs the AR 
nonlinear system’s output will remain bounded. This is 
because the order of the nonlinearity increases with time. 
[Ref. 12:pp. 68-69] 

Using this alternate tensor form of the AR nonlinear 
system, along with the generalized lattice of section 
II.D.4, the nonlinear lattice filter will be developed. To 
simplify the discussion, the filter’s finite memory will be 
restricted to N=2. Then the product Y(k) = cy eH1) 9% k-2) 1 
for A, ,A, =0,1,...,P forms a second order tensor, or ae two- 


dimensional data field given by: 
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E 


1 y(k-2) ees y (k-2) 
P 
y(k-1) y(k-1l)y(k-2) ... y(k-l)y (k-2) 
Yk =e 2 2 P (3 See 
y (k-1) y (k-1)y(k-2)..° y (h-))y eG 
P P P P 
y (k-1) y (k-l)y(k-2)... y (k-l1l)y (k-2) 


The types of systems that this nonlinear lattice 
structure can model exactly are of the form shown in Figure 
3.2 . The nonlinear combinations block forms a weighted sum 
of the cross-products and powers of the input variables 
y(k-i), for =o ae As an example, with N=2 and P=2 the 
general equation for the system’s output when excited by the 


input x(k) is given by: 


y(k) = x(k) - { H + H y(k-1) + H y(k-2) 


ile Zz 12 
Z Z 
+ H y(k-l)y(k-2) + H y (k-1) + HH y (k-1)y(k-2) 
a2 31 ee 
2 Z 
+ H y (k-2) + H y(k-l)y (k-2) 
ie 23 
C 2 
f° yo kaye VR (3 7s 


33 


It is also assumed that the system is time-invariant; 
otherwise the model changes with time k. In order to define 
a 

the forward and backward prediction errors, the (Pl) 


elements of the two-dimensional data field must by converted 


86 





XO 
YHNIGWOO 
UVANTINON 


seq 


(41) 


Figure 3.2 Nonlinear System Model 
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into a one-dimensional sequence. The ordering chosen by 


LENK is listed below: 


order = {(P,P),(P-1,P),(P,P=-1),..-,(0,2)3geee 
(P-1,P-1),(P-2,P-1),(P-1,P-2),...,(0,P-1), 


(P-1,0),...,(1,1), (0,1) 40,0) Gee 3 (3.15) 


where the components of the data field are identified by the 

indices (m,n)~ for m,n = 0,1, ee The elements of the 
2 

ordered set are numbered consecutively from 0 to (P+1) -1. 


The notation (m,n)-q is used to identify the q-th element 


prior to the element (m,n) as referenced to the ordered 


sequence. This notation will be further abbreviated to 

mn-q. Now the (q-1l)-order, normalized, forward error 

associated with predicting the value of the element 
m n 


y (k-l)y (k-2) from the preceding (q-1) elements is given by 


lee q-l A, A 
= a (myn)y (k-1)y (k-2) (3.16) 
mn 1A 
q-1 
for A, ,Az =0,1,2,...,P. Note that the coeffrterencs ay ,08n 
vAZ 
be thought of as the components of a second order’ tensor. 


q-1 
The coefficient ay 4 is equal to zero when the indices 
nova 


(A; Az) do not correspond to the (q-1) elements preceding 


(m,n) in the ordered sequence (i.e. when (A,,A3) >» (m,n) 
or when (A, ,A2) < mn-q ). Also, when (A\,Az) = (m,n), then 
q-l q-1l 
a (m,n) = 1 / Iti e Ly ae (3.17) 
mn mn 
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Similarly, the normalized backward prediction error in 
estimating y(mn-q) from the next (q-1) points is given by 
_q-1 q-t At Az 
r = 6 (mn-q)y (k-l)y (k-2) (3.18) 
mn-q ArAz 


q-l 
for Ai Ar = Oi, 1 ply. 5 Ps The coefficients Eten ad equal 
vAN2 


zero when the indices (A,,A,) do not correspond to the (q-1) 
elements following (m,n) (that is, when (),,A,) < (mn-q) or 
(Ay, A,) > (myn) ). When (A,,A2) = (m,n), then 
q-1l q-l1 

b (mn-q) = 1/ JIr lt. (3.19) 

mn-q mn-q 
(Ref. 12:pp. 145-146] 

Using the normalized nonlinear Levinson algorithn, 
LENK shows that the nonlinear prediction errors can be 


updated in order through the recursion relation 


_q _qa-1 
e e 
mn mn mn 
= §(K ) (3.20) 
_q q mle 
r r 
mn-q mn-q | 
where 
mn 
mn 1 1 - K 
OK )= —— q (ee) 
q mn 2 mn 
1 - (K ) -K 1 
q q 


89 


and the unique reflection coefficients are given by 


mn it ace 


K = Ef{e G ae (3.229 
q mn mn-q 


[Ref. 12:pp. 146-147] 


A deficiency remains to be corrected: the goal is 
to estimate y(k), but there is no y(k) term in the two- 
dimensional data matrix of equation (3.13). This problem is 


solved by adding another channel to the lattice structure 
and by exploiting the orthogonality of the backward pre- 
diction errors. Since the backward prediction errors leaving 
the top row of the lattice filter (Figure 2.9) are uncor- 


related, they can be used ina Fourier series to estimate 


y({k). These backward prediction errors can be formed into a 
2 
length L = (Ptl) vector defined as 
N 0 1 2 L-1 oT; 
[r i= oir i ae eo ] (3. 2e0) 
(m,n)-ZX 00 00-1 00-2 00-L+1 
where the subscript is in the form mn-q. Now, the error in 


estimating y(k) using the data in the Y(k) tensor is 


calculated from 


€= y(k)— kaa (3.24) 


where the Fourier coefficients, K.», are given by 


A 


20 


0 a | L-1 
[K,J=([E{y(k)r },E{y(k)r- prore, BLY CK) T }}.(3.25) 
r 00 00-1 OO = ter 
The resulting generalized nonlinear analysis lattice 


structure is depicted in Figure 3.3. [Ref. 12:pp. 147-150] 
The nonlinear lattice structure will now be examined 
more closely. As can be seen in Figure 3.3 , the inputs to 
the analysis filter are the normalized values of y(k) and 
the ry) 4 ordered components of the two-dimensional data 
field (equation (3.13)). The ordering of these inputs is in 
accordance with equation (3.15). At a given time k, these 
42 inputs are evaluated and inserted into the left 
side of the lattice structure. The lattice calculations are 
conducted by starting at the top row and moving downward 
along the northwest-southeast diagonal, and then advancing 
to the top of the next diagonal and so on. When all the 
lattice calculations for time k have been completed, the 
process is repeated for time k#tl. Just as in the case of 
the linear lattice filter, the PARCOR coefficients at each 
lattice section act to decorrelate the two input signals. 
The backward error signals exit at the top of the lattice 
Structure, and the forward error signals exit at the right 
Side. The output of interest for inverse filtering is the 
forward error signal from the top row of the filter (i.e., 


the row corresponding to the y(k) input signal). This is 


the error arising from the estimation of y(k). 
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Figure 3.3 Quadratic Nonlinear Lattice Filter 
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y¥ (k-1)y (k—2) 


Sr. The Nonlinear Lattice Applied to Deconvolution 


The generalized nonlinear lattice filter can be 
applied to the inverse filtering problem just as the linear 
lattice filter was in section II.D.6. The first problem is 
one of system identification and parameter estimation. Once 
the model parameters have been estimated, they are embedded 
in the nonlinear filter lattice structure. Then, the 
nonlinear lattice represents the inverse of the system’s 
transfer function. As will be shown, the nonlinear lattice 
filter can model both linear and nonlinear systems-- the 
linear system is just a special case of the nonlinear 
system. This inverse filtering algorithm was implemented by 
the FORTRAN programs NLMAIN, NLCLAT, SCHUR, NORMS, NLLAT, 
and URAND. These programs are listed in Appendix B. Now, 
this inverse filtering procedure will be described in more 
detail. 

In order to identify the model parameters (i.e., the 
PARCOR coefficients), the system was excited with zero mean, 
unit variance white noise sequences. Both Guassian and 
uniform noise distributions were used with good results. 
One difficulty encountered in generating the nonlinear data 
was ensuring that the output of the postulated nonlinear 
system remained bounded for the input noise signal. For the 
Simulations, the filter and system were constrained toa 


finite memory of at most two delays (N=2), while the 
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nonlinear order was allowed to vary from P=0’to P=4. (For N 
greater than two, the problem of ordering the elements of 
the N-dimensional Y tensor increases in complexity. ) The 
system output sequence y(k) was then processed by 
subroutines NLCLAT and SCHUR, which determined the 
corresponding autocorrelation matrix and the partial 
correlation coefficients, respectively. In an effort to 
improve the accuracy of these calculations, noise sequences. 
of up to 5,000 points were used. Additionally, the PARCOR 
coefficients were averaged over as many as 50 realizations 
of the input noise random process. Through trial and error, 
it was found that the best inverse filtering results were 
obtained when the resulting reflection coefficients were 
truncated to two decimal places. The output of subroutine 
SCHUR isa ets vs x Sie upper triangular matrix of 
reflection coefficients. This matrix 1s simply overlaid 
atop the upper triangular shaped nonlinear lattice structure 
to place the reflection coefficients at the correct filter 
sections. 

With the reflection coefficients calculated and 
embedded in the filter, we were able to use the lattice in 
an inverse filtering application. To recover the input 
Signal x(k), the system’s output signal y(k) was processed 
by subroutines NORMS and NLLAT. The function of NORMS was 


2 
to normalize the (P+1) +1 input signals into the lattice. 
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Then subroutine NLLAT, using the previously calculated 
reflection coefficients, implemented the lattice structure 
and carried out the lattice filter calculations. The forward 
error signal out of the top row of the lattice yielded the 
normalized estimate of x(k). Since the inputs to the lat- 
tice filter are all normalized, the outputs’ must be 
denormalized. Therefore, since the input into the top row of 
the lattice is divided by the norm of y(k), the forward 


error signal from the top row of the lattice is multiplied 


by the norm of y(k). (Note that if y(k) 1s a zero mean 
sequence, then this norm is equivalent to its standard 
deviation. ) Also, it was found that to achieve a good 
estimate, it was necessary to further scale this error 


Signal by dividing it by the norm of the noise generated 
sequence y(k). 

In order to verify the accuracy of the reflection 
coefficients, the noise generated output signal, y(k), was 
passed through the inverse lattice filter to see how well 
the filter whitened this sequence. The effectiveness of the 
whitening filter was evaluated by examining the mean-square 
error (MSE) between the white noise input signal, x(k), and 
the lattice filter’s output, Sau The running average 


mean-square error was calculated using the equation 





2 


kK A 
MSE(k) =\/(1/k) » (x (iy eo ey) (3.26) 
= | 





95 


Since the noise input, x(k), has unit variance, the MSE 
plot essentially provides a percentage error between x(k) 
and x(k). Pilots of the MSE are included in the 
simulation results. Having demonstrated that the lattice 
filter represented a good inverse of the system H(z), the 
system was then driven by a Known input signal x(k). To 
recover an approximation to this signal, the corresponding 
system output was passed eeraten the lattice filter, and . 
resulting lattice filter output was denormalized and 
rescaled as previously discussed. Simulation results for 
various systems ace presented in the aoilibergiun section. 
4. Inverse Filtering Simulation Results 

In this section, the previous modeling and inverse 
filtering procedures are implemented and applied to various 
linear and nonlinear systems. Unless otherwise noted, the 
reflection coefficients for each of the systems were 
determined by using twenty-five realizations (5,000 points 
each ) of the zero mean, unit variance white noise random 
process as the input excitation Signal. As previously 
mentioned, the mean-square error between this white noise 
input and the lattice filter’s output was plotted to 
evaluate the lattice filter’s inverse filtering performance. 
After the system was modeled, it was driven by the signal 
x(k) which consisted of ramps, pulses, and sinusoids as 


shown in Figure 3.4. 
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Figure 3.4 Input Signal x(k) 
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a. System I: y(k) = x(k) - (0.6y(k-1) +0.08y(k-2Z) ) 

This is the linear system first introduced in 

Section II.D.6. An input Gaussian white noise sequence was 
used to model the system. Figure 3.5 is a plot of. the 
running average mean-square error between the input noise 
and output error signals. In Section II.D.6, the linear 
lattice filter’s reflection coefficients were found to be 
K “=°-0.555 and K = '=0.077. It should be noted that these 
Sa appear in om first row of the nonlinear lattice’s 
reflection coefficient matrix. Here, for a first order, 


nonlinear lattice filter, the reflection coefficients are 


given by the upper triangular matrix 


0 QO -.55 -.07 0 
0 0 0 0 --.43 
K = 0 0 O -.55 0 
0 0 0 0 0 
0 0 0 0 0 


The system output, y(k), corresponding to the input signal 


of Figure 3.4 is shown in Figure 3.6. As can be seen by the 
plots of x(k) and % (ir) in Figure 3.7, the nonlinear lattice 
filter did an outstanding job of recovering the “unknown” 


input signal x(k) from the linear system’s output y(k). 
b. System II: y(k) = x(k) - (0.2y(k-1)y(k-2Z2)) 
This system involves a "cross-talk" 
nonlinearity, that is, the output is a function of the 


product of two different signals. For this system, both 
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Figure 3.5 System I Mean-Square Error 
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Figure 3.6 System I Output y(k) 
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Figure 3.7 System I Comparison of x(k) and x(k) 


Gaussian and uniform noise distributions were used in the 
modeling process in order to compare the two techniques. 


Here, both methods yielded the same reflection coefficients: 


ox 

tH 
Oo 2 2670 
oOoO0O0 © 
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oOo ore 
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The only difference between the two techniques was the value 


of the scaling factor (1.e., the norm of the noise generated 
System output y(k)). The plot of y(k) is shown in Figure 
SaOn The running average mean-square error and S(k) plots 


for the Gaussian noise derived model are depicted in Figures 
3.9 and 3.10, respectively. The corresponding plots for the 
uniform noise derived model are given in Figures 3.11 and 
3.12. As can be seen by comparing the plots, both modeling 
variations lead to nearly identical results. Both 
techniques yielded excellent approximations to the input 
Signal x(k). 
c. System III: y(k) = x(k) - O0.2y(k-1)y(K-1) 

This system involves a quadratic nonlinearity. 
Therefore, a second order model was used. The system output 
would not remain bounded for a Gaussian noise input, so the 
system was modeled using a unit variance uniform noise 


random process. The resulting reflection coefficients are: 
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Figure 3.9 System II Mean-Square Error 
(Gaussian Noise Model) 
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Figure 3.10 System II Comparison of x(k) and x(k) 
(Gaussian Noise Model) 
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(Uniform Noise Model) 
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Figure 3.12 System II Comparison of x(k) and S(k) 
(Uniform Noise Model) 
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Figure 3.13 provides the running average mean-Square error 
plot. The output y(k) of the nonlinear system is shown in 
Figure 3.14. Figure 3.15 depicts the ponponison between the 
System input x(k) and the lattice filter output K (ic). The 
x(k) curve has the same shape as x(k), however, it is 
slightly offset. In an attempt to improve the approximation 
process, the number of realizations of the noise random 
process used to model the system was doubled from twenty- 
five to fifty and the simulation was’ repeated. Although 
several reflection coefficient values changed by .01, there 
was no perceptible improvement in the estimation of x(k). 
d. System IV: y(k)=x(k) - (0.5+0.6y(k-1)+.08y(k-2) ) 
This example consists of the linear System I 
modified by the addidtion of a constant bias term. Gaussian 
noise was used to model the system. The reflection 


coefficients determined for the first order model are: 
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Figure 3.15 System III Comparison of x(k) and = (kK) 


0 (=623" 4 000 =. 0" 0 
0 O -.23 -.40 -.43 
K = 0 0 O =. 20 gee 
0 0 0 0 0 
0 0 0 0 0 


This is basically the same matrix as obtained for System I, 
but modified by the addition of several terms which attempt 
to model the constant offset. (Note the reappearance of the 
values of -.55 and -.07 in the topeere;) or Seen lattice 
matrix.) The running average mean-square error is shown in 
Figure 3.16. The system output y(k) is plotted in Figure 
Ss ealiegae Figure 3.18 is a plot of the lattice output Sk) 
and the system input x(k). The ee) curve 1S an excellent 
replica of the original input signal, but it is offset by a 
constant value. Although the small mean-Square error evident 
in Figure 3.16 indicates that the lattice is a good inverse 
filter, the lattice was unable to completely remove the bias 
term. Increasing the order of the model (1.e., overmodeling 
the system) did not improve the results. 
e. System V: y(k)=x(k) - (5.0+0.6y(k-1)+.08y(k-2) ) 
In order to further investigate the constant 
offset nonlinearity, the example of System IV was repeated 
uSing a bias of 5.0 instead of 0.5. Again, Guassian white 
noise was used in determining the model’s parameters. The 
reflection coefficients of the first order nonlinear lattice 


are given by: 
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Figure 3.16 System IV Mean-Square Error 


1s00 =6s- 2000s 2500 = 3000) 3500 39s 4000 =Sss 4500 = S000 
Pe earn le ere 


1000 


500 





O00S 


OOSF 


OO00F 


ACAIBMAIW TD PRIA AIDIAWIAONNDD 1 W OMADNMAONATASAS 


me ees — = 
: -—— « a 


baTWUS IW J 
O0Sd O00 


OO! 


00 


SZ 0- CO 


00°0 5C.0= a asGs 


CO 


BL e Wee 0S a0 


00° 


CM) A 
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System IV Comparison of x(k) and x(k) 
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0 =.92°=255 =2707 0 
0 O -.92 -.86 -.74 
Ka = 0 0 0 «(8B =D 
0 0 0 O -.90 
0 0 0 0 0 : 
Again, the reflection coefficient values of -.55 and -.07 


occur in the top row of the matrix as these elements 
apparently model the linear portion of the sytem. The 
running average mean-square enor is displayed in Figure 
3.19, and the output of the nonlinear system is shown in 
Figure 3.20. Figure 3.21 provides the comparison between the 
system input and the output of the inverse filter. As with 
System IV, the x(k) curve produced by the deconvolution 
process has the same shape as x(k), but is offset from the 
desired curve by a small constant value. Here, the system 
bias is 5.0, andthe x(k) and <(k) curves differ by about 
0.5. This is a relative improvement over the results 
obtained for System IV. System IV had a constant bias’ of 
only 0.5 but the x(k) curve was offset from the x(k) curve 


by approximately 0.4. 
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System V Output y(k) 


118 





9000 


4500 





S*] Olea S°0 00 6 Qe 
(4) LYHX“ (4X 


A 
Figure 3.21 System V Comparison of x(k) and x(k) 


119 


LEGEND 


NZ 
i 
Cs 
ale 
S| 


sre 


IV. CONCLUSION 


In this thesis, several linear least-squares 
deconvolution, Or inverse filtering, techniques were 
reviewed. Particular emphasis was placed on the lattice 
filter. It was shown that when a system is defined by an 
autoregressive, linear time-invariant model, this model can 
be transformed in to an equivalent lattice filter 
representation. Furthermore, the resulting analysis lattice 


filter acts as a whitening filter, and is the inverse of the 
system’s autoregressive transfer function. An example 
demonstrated the effectiveness of the lattice filter in a 
linear deconvolution application. Unfortunately, the linear 
lattice filter is unable to model nonlinear systems and, 
therefore, is not a viable nonlinear inverse filtering 
technique. 

The discussion of the linear lattice filter led to the 
development of the generalized lattice filter, which in turn 
led to the derivation of the nonlinear lattice filter. 
The goal of this thesis was to implement the nonlinear 


lattice filter, and then apply it to the linear and 


nonlinear deconvolution problem. This was accomplished with 
generally very good results. It was shown that the 
nonlinear batvtiuce filter was suitable for modeling 


ee 9 


discrete autoregressive systems where the output y(n) 
consisted of a weighted sum of the cross-products of 
the terms y (n-1) and Saino where the powers p and q can 
take on the values {0,1,2,3,4}. Examples were presented 
demonstrating the ability of the nonlinear lattice filter to 


effectively act as an inverse filter for both linear and 


nonlinear autoregressive systems. 
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APPENDIX A 
LINEAR LATTICE FILTER FORTRAN PROGRAMS 
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BEGIN CALCULATING THE REFLECTION FACTORS 


DO 50 J = 2,N 
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Noise 
Doel = 1 Ned 
he aie 
Dee 
RHO L,JI1) = ALPHA(I JT1)/BETA( IP JI } 
RNORM’ = DSQRT(I.0 - RHO(1I,JI1)*RHO(I,J11)) 
Doren kK = TON 
ie ALEHAC TSK 
ALPHAC T t = Bye ier Pala Fs anor ned NORM 
BETA(I,K) = (BETA(IP1,K)-RHO(I ,JI1)*T)/RNORM 
30 CONTINUE 
40 CONTINUE 
C WRITEC eae ees K=1,N),I=1,N) 
C42 FORMAT(/2X,13 (2x ; 2.5 ) 
i WRITE(8,42)J,((BETA(I,K) ,K=1,N),1=1,N) 
20 CONTINUE 
6 
RETURN 
END 
(AOI IS III IIIS OID OISISISISIICICIBIOIIOIOI SOIC ICICI ICICI IOI IOI IOI ICICI IOI III IK 
C 
C FUNCTION URAND 
C TAKEN _FROM COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS" BY 
C G.E. FORSYTHE, M.A. MALCOLM, AND C.B. MOLER 
TE cc ccnanxdkununuaniauavakanavthavaunwadiannexexnns 
C 
REAL FUNCTION LEO AEA 
INTEGER IA,IC,IT , MIC 
DOUBLE PRECISION H 
REAL S 
DOUBLE PRECISION DATAN,DSQRT 
DATA M2/0/ 1 1WO/3/ 
; IF(M2.NE.0) GO TO 20 
C IF FIRST ENTRY, COMPUTE MACHINE INTEGER WORD LENGTH 
M=1 
0 M2=M 
M=ITWO*M2 
TF(M. GT.M2) GO TO 10 
; HALFM=M2 
C COMPUT MULTIPLIER AND INCREMENT FOR LINEAR CONGRUENTIAL METHOD 
TA=8* IDINT( HALFM*DATAN( 1. D0)/8. pa} 
= 27 IDINT HALFM*(, 5D0-DSQR1T(3. D0)/6. DO) )+1 
; MIC=(M2-IC )+M2 
C S IS THE SCALE FACTOR FOR CONVERTING TO FLOATING POINT 
, S=. 5/HALFM 
C COMPUTE NEXT RANDOM NUMBER 
Zk AG Y= TY <1 
C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHICH DO NOT ALLOW 
C INTEGER OVERFLOW ON ADDITION 
: IF(IY.GT.MIC) IY= Si M2)-M2 
- IY=IY+IC 
C THE FOLLOWING IS FOR COMPUTERS FOR WHICH THE WORD LENGTH 
C FOR ADDITION IS GREATER THAN FOR MULTIPLICATION 
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IF(IY/2. GT. M2) TY=( 1Y-M2)-M2 
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LT ODI (1y+2)+M2 
FLOAT( IY) *S 
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TIME SEQUENCE IN AN ORDER 


WRITTEN 7 MAY 1985 BY P.J, LENK CRER atone) 


A 


HO 
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SUBROUTINE NLCLAT cy IS 
REAL*8 Y(5000) ,R(26,26) 


DEFINE CONSTANTS 
N* 


UN, ‘Yec(26) 


= N*N 
MNPL= MN + 1 
IYSM2 = [YS - : 

FIYSM2 = FLOAT(IYSM2) 
INITIALIZE R MATRIX TO ZERO 


oP 
es 
ae 


R(T.) = 0.0 
CONTINUE 
CONTINUE 
BEGIN OUTER LOOP 
DO 80 I = 3. 1YS 
becct ; Y(I) 
NSO MPL Sd N 
MO = MP1 21 
LIM = 2*MP1 - 1 
BO 40 b= 1 .LLIM 
iO = (eal 
Ti = MO 
31 = 10/2 
IF (MOD(LO,2).EQ.0) GO TO 30 
iS 
ji = MO 
IR TECCIR) = COORD(Y(I-1),I1)*COORD(Y(I-2) ,J1) 
CONTINUE 
CONTINUE 


CALCULATE THE CORRELATIONS 
DO 70 J = 1,MNPI1 
DO J 


a) Q) CA + a: (JYVEGCK) 
yRETECD 932 VEC9} ect 
CONTINUE 
CONTINUE 
CONTINUE 
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C DIVIDE BY THE NUMBER OF DATA ELEMENTS CONSIDERED 
DO 100 J = 1,MNP1 
00 ASK) 2 RUS K)/FIYSM2 
90 CONTINUE ’ 
100 CONTINUE 
C FILL IN THE SYMMETRIC HALF OF CORRELATION MATRIX 
DO 120 I = 2,MNP1 
ved 
HO eed =" om 
RCI) = 'R(J,1) 
110 CONTINUE 
120 CONTINUE 
C 
RETURN 
END 
C 
a KKK HH 4 x RRA A AAA RA REAR AR ARR RHI 
C 
C FUNCTION COORD 
C GENERATES OUTPUT OF RANDOM FUNCTION 
C CREATED 23 AUG 84 (REF 12:P. 187) 
FSC IONE UIEI ER ORIOLE IO RIE RIERA ERR RARER RRR 
DOUBLE PRECISION FUNCTION COORD/ x 1) 
C USE SIMPLE POWER SERIES TYPE POLYNOMIALS 
Y= 1.0 
IF (I. 9. 0) GO TO 30 
axe 
30 COORD = DBLE(Y) 
RETUR 
END 
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PREVIOUSLY CALCULATE 


WRITTEN 30 APRIL 86 
He He eK eK Ke eK KK kK Kk kk kK kK kk Kk kk kk ok ok kkk ok kk kk kk ok kok ok kk Kk 


SUBROUTINE Haig RHO ,N,NUMPTS ,NORM, XHAT ) 
KK IN R&R * 


EMENTS EXC E te haieR USING 
DeReeLe 


CITI MCI CCCI 9.190 CIO CI ONICOIC IC IAM C9 CI9C1IO8Y 


VARIABLE DEFINITIONS 
INFWD(1) = FORWARD ERROR INPUT INTO THE I'TH ROW OF THE LATTICE 
INBKD( I) = BACKWARD " ! 
OUT FW = FORWARD ERROR OUTPUT FROM THE I'TH ROW OF THE LATTICE 
OUTBKD I = BACKWARD rH i i J 11 I J J i 
Y = INPUT DATA VECTOR 
RHO = REFLECTION FACTOR MATRIX 
NORM = VECTOR OF NORMS OF THE LATTICE INPUT TERMS 
XHAT = OUTPUT DATA VECTOR; IT IS THE FORWARD ERROR SIGNAL FROM 
THE LAST STAGE OF THE FIRST ROW 

NUMPTS = NUMBER OF POINTS IN THE INPUT/OUTPUT SEQUENCES 
N = DIMENSION OF SQUARE Y DATA MATRIX 
MNP1 = DIMENSION OF THE RHO, NORM, AND INPUT/OUTPUT ARRAYS 
***VARTABLE DECLARATIONS**** 

INTEGER N,MN,NUMPTS, LAST ,MNP1 


(5000), INFWO(26) 
NORM 
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be Xt 


26 
BK 
TICE INPUTS FOR EACH Ear 
THIS TIME ITERATION 
R CASE WHEN K> 


COORD( Y(K-1) ,11)*COORD(Y(K-2) ,J1)/NORM( IR) 
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